# Code to run probit models and create figure 2 in EUP paper "Who is Powerful"
# Jonathan Slapin, July 2005

rm(list=ls(all=TRUE))  # clears everything in R's memory

library(foreign)       # library required to load stata dataset
library(MASS)


#  Load data
Amsterdam <- read.dta("C:/transtemp.dta")
attach(Amsterdam)

#  Run Probit for Model 1
model1 <- glm(am~b+dk+d+gr+e+f+irl+lux+i+nl+a+p+sf+s+uk+ep+com-1,family=binomial("probit"))
printmodel <- summary(model1)
print(printmodel)
# Set Expected Values for Model1 Denmark
z<- 5000 #number of sims

mu <- as.matrix(coefficients(model1))
sigma <- as.matrix(vcov(model1))
bsim <- mvrnorm(n = 5000, mu, sigma)

low <- z*.025
high <- z*.975
middle <-z*.5

x<- matrix(nrow=14,ncol=1)
k<- matrix(nrow=14,ncol=1)
for (i in 0:13) {
    j<-i+1
    x[j,]<-i*(1/7)-1
    }
for (i in 1:14){
    k[i,] <-i*(1/7)-1
    }
    
    
xaxis1 <- seq(1:14)
mid1 <-seq(1:14)
lci1  <-seq(1:14)
uci1  <-seq(1:14)

plot(xaxis1,mid1, type="n",
    xlab="number of member states in support coalition",ylab="power", ylim=range(-.5,.3))   

#Denmark
for (i in 1:14) {
            xsim1 <- as.matrix(c(x[i,],1,x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],0,0))
            xsim2 <- as.matrix(c(k[i,],-1,k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],0,0))   
            pred  <- pnorm(bsim%*%xsim1)
            predn <- pnorm(bsim%*%xsim2)
            fdif  <- pred-predn
            fdif <- sort(fdif)
            lci1[i] <- fdif[low]
            mid1[i] <- fdif[middle]
            uci1[i] <- fdif[high] 
            }
            lines(xaxis1,mid1,lty=1)
            

# Set values and graph UK


for (i in 1:14) {
            xsim1 <- as.matrix(c(x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],1,0,0))
            xsim2 <- as.matrix(c(k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],-1,0,0))   
            pred  <- pnorm(bsim%*%xsim1)
            predn <- pnorm(bsim%*%xsim2)
            fdif  <- pred-predn
            fdif <- sort(fdif)
            lci1[i] <- fdif[low]
            mid1[i] <- fdif[middle]
            uci1[i] <- fdif[high] 
            } 
            lines(xaxis1,mid1,lty=5)

# Set values and graph Italy


for (i in 1:14) {
            xsim1 <- as.matrix(c(x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],1,x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],0,0))
            xsim2 <- as.matrix(c(k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],-1,k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],0,0))   
            pred  <- pnorm(bsim%*%xsim1)
            predn <- pnorm(bsim%*%xsim2)
            fdif  <- pred-predn
            fdif <- sort(fdif)
            lci1[i] <- fdif[low]
            mid1[i] <- fdif[middle]
            uci1[i] <- fdif[high] 
            }
            lines(xaxis1,mid1,lty=6)
            
#Set values and graph Sweden


for (i in 1:14) {
            xsim1 <- as.matrix(c(x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],1,x[i,],0,0))
            xsim2 <- as.matrix(c(k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],-1,k[i,],0,0))   
            pred  <- pnorm(bsim%*%xsim1)
            predn <- pnorm(bsim%*%xsim2)
            fdif  <- pred-predn
            fdif <- sort(fdif)
            lci1[i] <- fdif[low]
            mid1[i] <- fdif[middle]
            uci1[i] <- fdif[high] 
            }
            lines(xaxis1,mid1,lty=2)

# Set values and graph France

for (i in 1:14) {
            xsim1 <- as.matrix(c(x[i,],x[i,],x[i,],x[i,],x[i,],1,x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],0,0))
            xsim2 <- as.matrix(c(k[i,],k[i,],k[i,],k[i,],k[i,],-1,k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],0,0))   
            pred  <- pnorm(bsim%*%xsim1)
            predn <- pnorm(bsim%*%xsim2)
            fdif  <- pred-predn
            fdif <- sort(fdif)
            lci1[i] <- fdif[low]
            mid1[i] <- fdif[middle]
            uci1[i] <- fdif[high] 
            }
            lines(xaxis1,mid1,lty=3)


# Set values and graph Germany

for (i in 1:14) {
            xsim1 <- as.matrix(c(x[i,],x[i,],1,x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],x[i,],0,0))
            xsim2 <- as.matrix(c(k[i,],k[i,],-1,k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],k[i,],0,0))   
            pred  <- pnorm(bsim%*%xsim1)
            predn <- pnorm(bsim%*%xsim2)
            fdif  <- pred-predn
            fdif <- sort(fdif)
            lci1[i] <- fdif[low]
            mid1[i] <- fdif[middle]
            uci1[i] <- fdif[high] 
            }
lines(xaxis1,mid1,lty=4)

legend=c("Denmark","Sweden","France", "Germany", "UK","Italy")
v=c(1,2,3,4,5,6)
legend(11,-.3,legend,lty=v)
